* Replication of Figure 8, panel a, using Golder marginal effects plot
* You will need to specify the directory in which the repliaction data are stored in the second line of code 

clear
cd "[Directory with replication data here]"

version 11.0
#delimit ;
log using "sim8a.log", replace;
set more off;

use "fordham_wp2025.dta", clear ;

logit multilat black time interactb ag interacta lynchings interactl cotton textapp, cluster(icpsr);

preserve;

set seed 339487731;

drawnorm SN_b1-SN_b10, n(10000) means(e(b)) cov(e(V)) clear;

postutil clear;
postfile mypost prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat diff_lo diff_hi 
            using "sim8a.dta", replace;
            
noisily display "start";

local a=0 ;
while `a' <= 18 { ;

{;

scalar h_time=8;
scalar h_cotton=0.26;
scalar h_textapp=4.7;
scalar h_black=20;
scalar h_ag=20;
scalar h_lynchings=30;
scalar h_constant=1;

    generate x_betahat0  = SN_b1*h_black 
                            + SN_b2*`a'   
                            + SN_b3*h_black*`a' 
                            + SN_b4*h_ag 
                            + SN_b5*h_ag*`a' 
                            + SN_b6*h_lynchings 
                            + SN_b7*h_lynchings*`a' 
                            + SN_b8*h_cotton 
                            + SN_b9*h_textapp
                            + SN_b10*h_constant;
                            
    generate x_betahat1  = SN_b1*h_black 
                            + SN_b2*`a'   
                            + SN_b3*h_black*`a' 
                            + SN_b4*(h_ag+15) 
                            + SN_b5*(h_ag+15)*`a' 
                            + SN_b6*h_lynchings 
                            + SN_b7*h_lynchings*`a' 
                            + SN_b8*h_cotton 
                            + SN_b9*h_textapp
                            + SN_b10*h_constant;
    
    gen prob0 =normal(x_betahat0);
    gen prob1=normal(x_betahat1);
    gen diff=prob1-prob0;
    
    egen probhat0 =mean(prob0);
    egen probhat1=mean(prob1);
    egen diffhat=mean(diff);
    
    tempname prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat diff_lo diff_hi;  

    _pctile prob0, p(2.5,97.5) ;
    scalar `lo0' = r(r1);
    scalar `hi0' = r(r2);  
    
    _pctile prob1, p(2.5,97.5);
    scalar `lo1'= r(r1);
    scalar `hi1'= r(r2);  
    
    _pctile diff, p(2.5,97.5);
    scalar `diff_lo'= r(r1);
    scalar `diff_hi'= r(r2);  
   
    scalar `prob_hat0'=probhat0;
    scalar `prob_hat1'=probhat1;
    scalar `diff_hat'=diffhat;
    
    post mypost (`prob_hat0') (`lo0') (`hi0') (`prob_hat1') (`lo1') (`hi1') 
                (`diff_hat') (`diff_lo') (`diff_hi');
   
    };      
    
    drop x_betahat0 x_betahat1 prob0 prob1 diff probhat0 probhat1 diffhat;
    
    local a=`a'+ 1;
    
    display "." _c;
    
} ;

display "";

postclose mypost;

restore;

merge using "sim8a.dta";

gen yline=0;

gen MV = _n-1;

replace  MV=. if _n>18;

gen MV2=MV+1945;

graph twoway hist year, percent color(gs14) yaxis(2)
        ||  line diff_hat MV2, clwidth(medium) clcolor(blue) clcolor(black)
        ||  line diff_lo MV2, clpattern(dash) clwidth(thin) clcolor(black)
        ||  line diff_hi MV2, clpattern(dash) clwidth(thin) clcolor(black)
        ||  line yline MV2,  clwidth(thin) clcolor(black) clpattern(solid)
        ||  ,   
            xlabel(1945 1950 1955 1960 1962, nogrid labsize(3)) 
            ylabel(-0.2 -0.1 0 .1 .2, axis(1) nogrid labsize(3))
            ylabel(0 5 10 15 20 25, axis(2) nogrid labsize(3))
            yscale(noline alt) yscale(noline alt axis(2)) xscale(noline) legend(off)
            yline(0, lcolor(black))
            yline(-0.2 -0.1 0 .1 .2, lcolor(white))
            xtitle("Year", size(3))
            ytitle("Marginal Effect of 1 SD Increase in Agricultural Employment", size(3))
            ytitle("Percent of Observations" , axis(2) size(3))
            xsca(titlegap(4)) ysca(titlegap(4)) ysca(axis(2) titlegap(4))
            scheme(s2mono) graphregion(fcolor(white) ilcolor(white) lcolor(white));
            
log close;
